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ABSTRACT 

In the past year, the HiRes and Auger collaborations have reported the discovery of 
a high-energy cutoff in the ultra-high energy cosmic-ray (UHECR) spectrum, and an 
apparent clustering of the highest energy events towards nearby active galactic nuclei 
(AGNs). Consensus is building that such ~ lO^^-lO^*' eV particles are accelerated 
within the radio-bright lobes of these sources, but it is not yet clear how this actually 
happens. In this paper, we report (to our knowledge) the first treatment of stochastic 
particle acceleration in such environments from first principles, showing that energies 
~ 10^" eV are reached in ~ 10^ years for protons. However, our findings reopen the 
question regarding whether the high-energy cutoff is due solely to propagation effects, 
or whether it does in fact represent the maximum energy permitted by the acceleration 
process itself. 

Key words: cosmic rays - physical data and processes: acceleration of particles; 
plasmas; turbulence ~ galaxies: active; nuclei 



1 INTRODUCTION 

Cosmic rays are energetic charged particles traveling 
throughout the Galaxy and the intergalactic medium under 
the influence of various physical processes, including deflec- 
tion by magnetic fields, and collisions with other particles 
along their trajectory. Their energy spectrum measured at 
Earth is a steep (roughly power-law) distribution with log- 
arithmic index a ~ 2.6-3, extending up to a few times 10^° 
eV. Other than their spectrum, these particles are charac- 
terized by their angular distribution in the sky, and by their 
mass composition. 

A highly significant steepening in the UHECR spec- 
tru m was reported by both the HiRes (Abbasi et al. 2008) 
and lAugeil (|2008 j ) collaborations. (Distinguished from their 
lower-energy counterparts, UHECRs have energies in excess 
of 1 EeV = 10^® eV.) This result may be a strong con- 
firmation of the predicted Greisen-Zatsepin-Kuzmin (GZK) 
cutoff due to photohadronic interactions between the UHE- 
CRs and low-energy photons in the cosmic microwave back- 
ground (CMB) radiation (Greisen 1966; Zatsepin & Kuz'min 
1966). Together with the measured low fraction of high- 
energy photons in the CR distribution, this measurement 
already rules out so-called top-down models, in which the 
UHECRs represent the decay products of high-mass dark 
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matter particles created in the early Universe (Semikoz et 
al. 2007). The measured photon fiux is also in confiict with 
scenarios in which UHECRs are produced by collisions be- 
tween cosmic strings or topological defects (Bluemer et al. 
2008, Auger 2008b). On the other hand, such energetic par- 
ticles may still be produced via astrophysical acceleration 
mechanisms (see Torres & Anchordoqui 2004 and other ref- 
erences cited therein). 

UHECRs are not detected directly, but through the 
showers they create in Earth's atmosphere (see, e.g., Melia 
2009). Depending on the energy and type of primary par- 
ticle, the ensuing cascade has characteristics that allow the 
ground-based observatories to determine not only whether 
the incoming UHECR is a photon, but also its atomic num- 
ber. It should be pointed out, however, that a determination 
of the primaries' composition strongly relies on an extrapola- 
tion of current phenomenological hadron ic interaction mod- 
els, so it remains rather uncertain. The lAugeil l|2007l ) data 
confirmed the dominance of protons in primary cosmic rays, 
though they also exhibit evidence for a mixed composition 
extending to energies as high as ~ 50-60 EeV, with a higher 
atomic number Z, up to Z ~ 26 (Unger et al. 2007). 

But the most telling indicator for the possible origin of 
these UHECRs is the discovery by Auger (Auger 2008a) of 
their clustering towards nearby (~ 75 Mpc) AGNs along the 
supergalactic plane. The significance of this correlation has 
been further strenghtened by a more recent analysis which 
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weights the AGN spatial distribution by their hard X-ray 
flux (George et al. 2008) . This raises at least two questions: 
(1) How are the UHECRs accelerated to such high energies? 
and (2) given these nearby sources, is the sharp suppression 
of UHECRs in the last decade of their observed energies 
really due solely to the GZK effect, or does it signal a limi- 
tation to the acceleration efficiency? 

Previous attempts at understanding how particles are 
accelerated to EeV energies and beyond have generally been 
based on first-order Fermi acceleration (see, e.g., Ostrowski 
2008, and other references therein) within shocks created by 
blast waves like those in supernova remnants (Fatuzzo & 
Melia 2003, Crocker et al. 2005). But this process is sub- 
ject to kinematic restrictions that inhibit the particles from 
actually reaching ultra-high energies (see, e.g., Nayakshin 
& Melia 1998, Gallant et al. 1999). Recent numerical sim- 
ulations have shown that an increase in the Lorentz fac- 
tor 7 of ultra-relativistic shock waves steepens the observed 
spectrum (Niemiec & Ostrowski 2006) and reduces its high- 
energy cutoff. 

For these reasons, it is not plausible for UHECRs to 
emerge from astrophysical environments, such as supernova 
remnants, where first-order processes are dominant so long 
as the shock velocity is super-Alfvenic, because they cannot 
even contain such high-energy particles (Hillas 1984) — the 
gyration radius of particles with energy ~ 10^'' eV for a 
typical galactic magnetic field is much larger than the size 
(< 10 pc) of these structures. 

On the other hand, a second-order Fermi process (Fermi 
1949) can explain observational features not addressed by 
the first-order process, as in the case of a supernova rem- 
nant itself (see Cowsik & Sarkar 1984). Moreover, stochas- 
tic particle acceleration through a gyroresonant interaction 
with MHD turbulence (a second-order Fermi process; see 
Fermi 1949) can be very efficient if the Alfven velocity ap- 
proaches c (Dermer & Humi 2001). The stochastic accel- 
eration of particles by turbulent plasma waves has already 
received some attention in the literature (see Liu et al. 2004, 
2006, and references cited therein, and Wolfe & Melia 2006). 
Indeed, the feasibility of second-order Fermi acceleration in 
radio galaxies has been demonstrated through the steady 
re-acceleration of electrons in certain hot spots (Almudena 
Prieto et al. 2002). 

Our treatment from first principles, however, avoids 
many of the previously encountered unknowns and limi- 
tations. In this paper, we report (to our knowledge) the 
first treatment of stochastic acceleration of charged parti- 
cles in the lobes of radio-bright AGNs by directly comput- 
ing the trajectory of individual particles. An earlier version 
of this treatment — for the propagation of charged particles 
assumed already accelerated at TeV energy through the tur- 
bulent magnetic field at the Galactic centre — may be found 
in Ballantyne et al. 2007, and Wommer, Melia & Fatuzzo 
2008; by contrast, in the present paper both the propaga- 
tion and acceleration are taken into account. We show that 
random scatterings (a second-order Fermi process) between 
the charges and fiuctuations in a turbulent magnetic field 
can accelerate these particles up to ultra-high energies, pro- 
vided a broad range of fluctuations is present in the system. 



2 DESCRIPTION OF THE MODEL 

In our treatment, we follow the three-dimensional motion 
of individual particles within a time-varying turbulent mag- 
netic field. By avoiding the use of equations describing sta- 
tistical averages of the particle distribution, we mitigate our 
dependence on unknown factors, such as the diffusion co- 
efficient. We also avoid such limitations as the Parker ap- 
proximation (Padmanhaban 2001) in the transport equa- 
tion. However, a remaining unknown is the partitioning be- 
tween turbulent and background fields. For simplicity, we 
take the minimalist approach and assume that the magnetic 
energy is divided equally between the two components. 

Another unknown is the turbulent distribution. For 
many real astrophysical plasmas, the magnetic turbulence 
seems to be in accordance with the Kolmogorov spectrum. 
This is seen, e.g., in the solar wind (Leamon et al. 1998) 
and through interstellar scintillation (Lee & Jokipii 1976); a 
more recent numerical analysis of MHD turbulence confirms 
the general validity of the Kolmogorov power spectrum (Cho 
et al. 2003). In addition, renormalization group techniques 
applied to the analysis of MHD turbulence also favour a 
Kolmogorov power spectrum (for more details, see Smith et 
al. 1998, and Verma 2004). 

We model the radio lobe of an AGN as a sphere of radius 
TZ, a second parameter in our simulations. A population of 
relativistic particles of mass m, protons or heavy ions, with 
an energy E = ^mc? , where 7 is the Lorentz factor, is re- 
leased in an inner sphere of radius TZ' ~ aTZ. The value of a 
must be much smaller than 1, otherwise very few particles 
reach an energy E > lO'^^eV. For a small value of a, the 
gyration radius becomes comparable to the size of the accel- 
eration region at _E > 10^* eV, and therefore changes in a do 
not significantly alter the result. For the sake of specificity, 
we use a value a ~ 10"^ in this paper. Once released, the 
particles propagate through the turbulent field until they 
escape the acceleration region and enter intergalactic space. 

2.1 Time varying turbulent field 

We use the lGiacalone k, Jokipiil (Il994l ') prescription for gene- 
rating the turbulent magnetic field. Their principal aim of 
propagating individual particles through a magnetostatic 
field was to compute the Fokker-Planck coefficients for a 
direct comparison with analytic theory. For our purpose, we 
modify that prescription to include a time-dependent phase 
factor that allows for temporal variations. 

The global magnetic field is written as a sum of a back- 
ground term Bo, constant and uniform, and a turbulent field 
varying in space and time (i.e., as a superposition of Alfven 
waves) . 

The equation of motion of a relativistic test particle 
with charge e and mass m moving in an electromagnetic 
field F'^^ is the Lorentz equation (Landau & Lifchitz 1975, 
Meha 2001) 

du^ e 

mc = -F'^'^Ui, (1) 

ds c 

(with = 0,1,2,3), where c is the speed of light in vac- 
uum, u'^ — {"f,-yv/c) is the four-velocity of the particle, 
7 = 1/-^1 — {v/cY is the Lorentz factor, and s/c is the 
proper time. We calculate the trajectory of the particle in a 
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magnetic field B(t, r) = {mc/e)fl{t,r) as a solution of the 
space components {n — 1, 2, 3) of Equation (1) 



(2) 



where t is the time in the rest frame of the acceleration 
region. The quantity fl{t,r) in Equation ([2} is given by 

n{t,r) = no + sn{t,r) , (3) 

where flo = (e/mc)Bo, in terms of the background mag- 
netic field Bo, and Sfl{t, r) is the time-dependent turbulent 
magnetic field. We ignore any large-scale background electric 
fields — a reasonable assumption given that currents would 
quench any such fields within the radio lobes of AGNs. The 
time variation of the magnetic field, however, induces an 
electric field S£{t, r) = (e/mc)E(t, r) according to Faraday's 
law. 

The procedure of building the turbulence calls for the 
random generation of a given number N of transverse waves 
ki, i — 1,..,N at every point of physical space where the 
particle is found, each with a random amplitude, phase and 
orientation defined by angles 9{ki) and (jiiki). This form of 
the fiuctuation satisfies V • 5Q,{t,r) = 0. We write 



where the polarization vector is given by 
^±(ki) = cos a{ki)y' ± i sin a{ki)z . 



(4) 



(5) 



Given the form in Equation (4) for the turbulence, the elec- 
tric field 5£{t,r) is given by 



5£(t,Y) 



with 

^±(fci) = ±isina(fci)y' — cosQ(fci)z' 



(6) 



(7) 



The orthonormal primed coordinates r' = {x',y',z') are re- 
lated to the lab-frame coordinates r = {x, y, z) via the ro- 
tation matrix R{9,cj)), in such a way that for every k the 
propagation vector is parallel to the x' axis. The matrix 
R{0, (j}) is given by 



/ cos o cos (p 
r' = I — sin 4> 
\ — sin 6 cos ( 



cos u sm (p 

cos (j) 
— sin 6 sin ( 




(8) 



For each value of ki, there are 5 random numbers: < 
e{ki) < TT, < (l>iki) < 27r, < a{ki) < 2-k, < f3{ki) < 2-k 
and the sign plus or minus indicating the sense of polariza- 
tion. 

Further assumptions are necessary to specify the dis- 
persion relation uj — uj{ki). For every turbulent mode, we 
use the dispersion relation for transverse non-relativistic 
Alfven waves (see Kaplan & Tsytovich 1973 for an extended 
discussion): u!{ki) — VAkiCos6{ki), for i — 1,..,N, where 
va = -Bo / \/47i^'Ti^ is the non-relativistic Alfven velocity 
in a medium with background magnetic field Bo and num- 
ber density n, nip the proton mass, and 9{ki) is the angle 
between the wavevector ki and Bo. This is the condition 
thought to be valid for the propagation of turbulent modes 



in a magnetized astrophysical environment, such as the ra- 
dio lobes of an AGN. The background plasma is assumed to 
have a background proton number density n ~ 10~* cm~^, 
a reasonable value for these environments (Almudena Prieto 
et al. 2002). 

The amplitudes of the magnetic fiuctuations are as- 
sumed to be consistent with Kolmogorov turbulence, so 



^(^0 — ^{kmin^ { -J 1 



-r/2 



(9) 



for i — 1, .., A*', where kmin corresponds to the longest wave- 
length of the fluctuations and the index F of the power spec- 
trum f2^(fc) is 5/3. Finally, the quantity Q,{kmin) is com- 
puted by requiring that the energy density of the magnetic 
fiuctuations equals that of the background magnetic field: 



S = 



B^{k^ 



i=l 



Stt 



2 2 ^ 

^-n^(k )^ 



87re2 



J Stt 



(10) 



We choose N=2400 values of k evenly spaced on a logarith- 
mic scale; i.e., a wavenumber shell with bounds ki — ki+i 
holds fci+i — ki X {kmax/kmin)^^'^ values. Considering that 
the turbulence wavenumber k is related to the turbulent 
length scale I by k = 2n/l, we adopt a range of lengthscales 



from Imin = 10 ^ vo /Qo to Z„ 



10 Vo/Ho, where vo 



the initial velocity of the particle and Qo is the initial gy- 
rofrequency in the background magnetic field. Thus the dy- 
namic range covered by k is kmax/kmin — Imax/lmin = 10^°, 
and our description allows for 240 transverse modes k per 
decade. The values of kmax and kmin fix the magnetic en- 
ergy equipartition through Equation (10). The value km.in~^ 
is proportional through a factor of order 1 to the correlation 
length of the turbulence (see Ruffert and Melia 1994, and 
Rockefeller et al. 2004, for examples of how this is generated 
in the interstellar medium); the value kmax~^ is the wave- 
length at which the interaction between the turbulence and 
most of the particles is the most efficient, so that energy is 
drained out of 5Q. 

Ifowever, since the gyroradius rg{E) evolves over a 
large energy interval, the gyroresonant wavenumber kres (E) 
moves accordingly in such a way that in the global wavenum- 
ber interval {kmin ~ kmax) there is for every E a certain 
kres{E) fulfilling the resonance condition rg{E)kres{E) ~ 1. 
Such a k range involves a large computational time, espe- 
cially if a statistically significant number of particles is to 
be considered. 

Since our numerical simulation is not performed by 
specifying the magnetic turbulence on a computational grid 
with given cell size Ax, the choice of kmax = 2n/l,nin is not 
dictated by a fixed spatial resolution (see section |3] for more 
details). In addition, the result is not affected by spurious 
effects to the discreteness or the periodicity. As a bypro- 
duct, the divergenceless condition V • 5Cl{t,r) = is easily 
satisfied and does not require an extension of the Godunov 
solver of the MHD equations for the purpose of "divergence 
cleaning" (Ryu et al. 1998) or a reformulation of the MHD 
equations including, e.g., divergence-damping terms (Ded- 
ner et al. 2002). 

With this prescription, we construct the turbulent mag- 
netic field at every point of physical space where the particle 
is found, which we then propagate without taking any time- 
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average along the trajectory. The particles passing through 
this region are released initially at a random position inside 
the acceleration zone, which for simplicity is taken to be a 
sphere of radius TZ, with a fixed initial velocity vo pointed 
in a random direction. T he initial value of the Lorentz fac- 
tor 70 — — {vq/cY — 1.015 is chosen to avoid having 
to deal with ionization losses for the protons or heavy ions. 
The particles closest to the edge of the acceleration zone 
have a higher probability of escaping than those starting 
farther in, and therefore reach relatively lower energies. In 
the usual (Fermi) way, this produces (in the highest energy 
portion of the spectrum) an inverted power-law distribution. 



2.2 Energy losses 

In principle, energy losses due to synchrotron and inverse 
Compton processes involving radio and Cosmic Microwave 
Background (CMB) photons, all of which increase as 7^, 
can significantly limit the maximum energy attainable by 
a cosmic ray during the acceleration process, given that its 
Lorentz factor 7 evolves from ~ 1 up to 10^" — 10^^. For an 
UHE particle (either a proton or a heavy ion) , both the radio 
and CMB photons will have an energy 'yhv in the centre-of- 
momentum frame well below the rest energy of the cosmic 
ray (i.e., ^hv « mc^ , where m is the mass of the accele- 
rating particle). For the purpose of these estimates, we use 
a radio frequency Vradio = 0.1 GHz {huradio ~ 4.2 x 10~^ 
eV) and a CMB frequency corresponding to the peak of 
the blackbody spectrum, vcmb = 2.821kT/h — 158 GHz 
{hvcMB ~ 6.6 X 10~* eV), where k is the Boltzmann con- 
stant, h the Planck constant, and T = 2.7 K is the CMB 
temperature. Consequently, the energy losses due to inverse 
Compton may be calculated in the Thomson limit. Com- 
pare this with the situation for high energy electrons, for 
which the Thomson condition would not be satisfied even 
at energies ^m^cp' ~ 10^^ eV, requiring in that case the full 
Klein-Nishina treatment. 

The propagation of high-energy particles is here mo- 
deled in a region of tens of kpc size. Therefore we neglect 
any effect of the relativistically-narrowed jet on the spatial 
distribution of the radio background, assumed for simplic- 
ity to be isotropic. Since the CMB intensity field is also 
isotropic, we take these energy losses into account using the 
following angle-integrated power-loss rate: 

- = jc^T(m)c'y^ I— + Ur + Ucmb] , (11) 

where ctt (m) = 6.6524x(me/m)^ lO"^'' cm'^ is the Thomson 
cross section for a generic particle of mass m, which can 
be a proton or heavy ion, and i3^/(87r) = (2_Bo'^)/(87r) is 
the total energy density of the magnetic field. The photon 
energy density Ur inside a typical radio lobe is computed 
as Ur = L/{4-KcR.^), where we assume L to be a standard 
luminosity density corresponding to the Fanaroff- Riley class 
II of galaxies (i' = 5 x lO^'^ W Hz~^ sr^^ at 178 MHz), 
and TZ is the size of the spherical acceleration zone. For the 
CMB, we use Ucmb = aT'^ = 4.2 x 10"" erg cm"^. In 
a region where magnetic turbulence is absent or static, a 
given test particle propagates by "bouncing" randomly off 
the inhomogeneities in B, but its energy remains constant. 
The field we will model below, however, is comprised of time- 




Figure 1. Three- dimensional trajector y of sin gle particle in the 
turbulent field of iGiacalone fc Jokipiil l ll994l ) reproduced with 
our code. The length scales are in units of the gyration radius, 
i>o/f^O) which remains constant during the propagation. The en- 
ergy is verified to be constant, as expected, over a time inter- 
val At = lOOOOo"^, within a relative error of 10~^. The back- 
gr ound magnetic field Bp is p arallel to the z axis and, as found 
bv*Giacal one fc Jokipiil l ll994l) . the difi'usion along Bq dominates 
with respect to that across the field. 




Figure 2. Three-dimensional trajectory of a single particle re- 
leased at random within the acceleration zone, assumed to be 
a sphere of radius 50 kpc. The scale on the axes is in units of 
W^vo/Qq, where vg/i^o is the initial gyroradius. The particle is 
released with a fixed initial speed vq, but pointed in a randomly 
chosen direction. The calculation stops when the particle leaves 
the radio lobe and is injected into the intergalactic medium. 



varying gyroresonant turbulent waves (see Equation |4]), and 
collisions between the test particle and these waves produces 
a net acceleration (in the lab frame). 



3 NUMERICAL CODE SETUP 

In this section we describe the numerical code used to per- 
form the simulations. The Lorentz Equation ((2)1 is integrated 
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using a Runge-Kutta 4"* order method (Press at al. 1997) 
for the system of 6 first-order differential equations 



dxi 
"dT 
dui 



c 

-Ui 

7 

S£, - 



[u X n], 

7 



(12) 
(13) 



for i = 1,2,3. The components S£i and ili are intended to 
be the real parts of the corresponding complex quantities. 

The portable random number generator used to produce 
the turbulence is Knuth's subtractive routine ran3 (Press et 
al. 1997), with a seed number / = lO". This routine has 
a relatively short execution time and is suitable to avoid 
the introduction of unwanted correlations into the numerical 
computation. 

There are two approaches to numerically implementing 
a turbulent magnetic field generated by this method. The 
first approach (used in Giacalone & Jokipii 1994) is to cal- 
culate the magnetic field at every time step for each particle 
position. The position is then found by solving the Lorentz 
Equation ([2]). In the second approach, the magnetic field is 
generated for a given volume at the beginning of the simu- 
lation and then it evolves according to Equation @. In or- 
der to have an acceptable fc-binning with a dynamical range 
of kmax/kmin = 10^°, oue would then need to specify the 
field at an excessively large number of lattice points. This is 
not only time-consuming, but also very memory-intensive. 
So, like Giacalone & Jokipii (1994), we adopt the former 
approach. In this way, the magnetic field is generated only 
where needed, and the overwhelming amount of computer 
memory required by the second approach is not necessary,. 
Since the confinement volume is a parameter of the model, 
the second approach would also require adapting the lattice 
spacing in order to maintain the same space resolution in 
physical space. 

The Runge-Kutta integrator has previously been vali- 
dated for the cases of uniform and constant electric and 
magnetic fields, where the outcome of the simulation can be 
compared with an analytical solution. 

Secondly, as a validation test of the code, we reproduced 
the result of Giacalone & Jokipii (1994) for the case of a 3D 
magnetostatic turbulence, by using the same set of para- 
meters. We discretize the turbulence in 50 transverse modes 
k, where the values of k are chosen to be evenly spaced in 
logarithmic scale in the interval of the corresponding length 
scales from Imin ~ (l/5)'i;o/no to Imax = lOvo/rio, where 
vo is the initial velocity of the particle and Qq is its gy- 
rofrequency in the background magnetic field. The particles 
are released in a random initial position with initial velocity 
randomly oriented but with fixed value vq. In Figure [T] we 
present the trajectory of a single particle, the position along 
the three axes expressed in units of the gyration radius. In 
this test, the energy of the particle is constant within a rel- 
ative error of 10~^ over a time interval corresponding to 

io^n-\ 

In order to produce the time-dependent turbulent mag- 
netic field considered in this paper, between two successive 
shufflings of all five random quantities in 5 ft, which are 
performed every At ~ 10* — lO^s, the particle propagates 
gyroresonantly with the oscillating turbulence. We verified 
that a change in the Runge-Kutta time-binning by over one 
order of magnitude does not produce any systematic numer- 



J lo" 




Time (s) 



Figure 3. Simulated time evolution of the Lorentz factor 7 for 
a proton propagating through a time- varying turbulent magnetic 
field. The three curves correspond to three different values of Bq: 
10"^, 10"®, and lO"'* gauss. The protons are released at an initial 
random position inside the acceleration zone — a sphere of radius 
7^ = 50 kpc — with the same initial speed vq, though pointed 
in random directions. The proton is followed until it leaves the 
acceleration zone and enters the intergalactic medium. The acce- 
leration timescale At is inversely proportional to the background 
field Bq. Therefore, as expected, a larger Bq produces a more 
efficient acceleration. In this example, a proton winding its way 
through a field Bq = 10~* gauss attains an energy E ~ 10^" eV 



in approximately 10 years. 



ical effects associated with y{t) and the spectrum. Changing 
the fc-binning from N = 1200 to iV = 3000 in Equation g] 
similarly does not noticeably change the resultant 7(t) and 
the spectrum. We chose A'^ — 2400 which results in a reason- 
ably long computational time. However, a coarser fc-birming, 
e.g., with 10 modes/decade, could possibly result in a worse 
determination of the macroscopic indicator as instantaneous 
spatial diffusion coefficients; the time evolution of the diffu- 
sion coefficients is however beyond the scope of the present 
paper. 



4 RESULTS AND DISCUSSION 

Figure [5] (to be compared with Figure [TJ shows the tra- 
jectory of a single particle released at random within the 
acceleration zone, with the initial speed vq, pointed in ran- 
dom direction, with an ambient magnetic field Bo = 10~* 
gauss. In Figure |3] we plot the time evolution of the par- 
ticle Lorentz factor 7 for three representative values of the 
background field Bo; 10^^, 10~*, and 10~^ gauss. We see 
the particle undergoing various phases of acceleration and 
deceleration as it encounters fluctuations in B. 

The acceleration of the particle results from the 0-th 
component of Equation [T] which reads 



7 



(14) 



where SE^Ui is the scalar product of the electric fleld and the 
3- velocity of the particle. Therefore the acceleration is given 

by 7(i) = \/ 70^ + 2 SS^Uidt'. The integrand function can 
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Energy (eV) 



Figure 4. Calculated differential spectra for 500 protons with 
Bq = 10~* gauss and for different values of the size of the 
acceleration region, assumed to be a sphere of radius spanning 
the interval = [5 — 160] kpc. The dependence of the energy 
cutoff on TZ is evidenced. This result shows that the cutoff in 
the observed spectral distribution can be due to the competi- 
tion between two distint effects: propagation through the CMB 
and intrinsic properties of the accelerator. Moreover, the slope 
in the region E > 4 X 10^^ eV strongly depends on R. This di- 
agram supports the view that the steeper CR spectrum below 
log{E /eV) ~ 18.6 likely represents a population of galactic cos- 
mic rays. 



be strongly time- varying and therefore needs to be computed 
numerically. 

As can be seen in Figure [S] the acceleration timescale 
Atacc is inversely proportional to Bo so, as expected, more 
energetic turbulence accelerates the particles more effi- 
ciently, in agreement with what was expected from the non- 
relativistic Alfven wave theory. Previous studies (Casse et 
al. 2002) of particle transport through a turbulent magnetic 
field using the prescription in Giacalone & Jokipii (1994), 
compared with a Fast Fourier Transform method, showed 
that the time of confinement within the jets of FR II galax- 
ies is too small for the particles to attain an energy of 10^" 
eV. In our case, the particle acceleration takes place over 
a much bigger volume (with dimension 7?. = 50 kpc, com- 
patible with the known size of FR-II radio galaxies) and 
the efficiency is enhanced by the strong temporal variation 
of the turbulence. We find, in particular, that particles can 
easily accelerate to UHE on timescales short compared to 
the age of the radio-lobe structure through a gyroresonant 
interaction with a magnetic turbulence. 

In our simulation, the particle acceleration is efficient 
because it occurs over a wide range of turbulent fluctua- 
tions, such that the wave-particle interaction is resonant at 
all times. Such a distribution is expected if the magnetic en- 
ergy cascade proceeds (without loss) from the largest spa- 
tial scales down to the region where energy dissipation and 
transfer to the particles becomes most efficient. 

In Figures 3] and \5\ we show the spectral distributions 
for distinct values of the turbulent energy and size of the 
acceleration region. In order to produce this result, we fol- 
lowed the trajectory of 500 protons, launched in the manner 
described above, with different values of the parameters R 
and Bq. 




1.0x10^^ 



Energy (eV) 



Figure 5. Calculated differential spectra for 500 protons with 
-R = 50 kpc and for different values of the turbulent magnetic 
energy. In this case Bq spans the interval Bg = 1.5 X [10~^ — 10~*] 
gauss. See the comments in Figure |4] 



Bg=6 X 1 0" gauss - 
power law with index -2.5 - 




5.0x10 1.0x10 1.5x10 2.0x10 2.5x10 3.0x10 

Energy (eV) 

Figure 6. Calculated differential spectrum for 1,000 protons in 
the energy range log(i?/eV') = [18.6—19.5] for the selected param- 
eters Bo = 10~* gauss and 7J. = 50 kpc. A power-law behaviour 
with index —2.6 in the differential spectrum of protons injected 
into the intergalactic medium in this model is in agreement with 
a recent statistical analysis of HiRes observations (Gelmini et al. 
2007). This good match supports the view that the steeper CR 
spectrum below log(i?/ey) 18.6 represents a different popula- 
tion, possibly associated with the Galaxy itself. 



We conclude that the observed spectral cutoff (Abbasi 
et al. 2008, Auger 2008b) can result from the competition of 
two distint effects: not only the GZK cutoff, namely degra- 
dation of primary UHECRs due to the propagation through 
the CMB, but also, and possibly dominant, intrinsic proper- 
ties of the source which constrain the process of acceleration. 

Figure |6] depicts the calculated differential spectrum in 
the energy range log(i?/eV) = [18.6 — 19.5]. From our sam- 
pling of the various physical parameters, we infer that for a 
radius TZ = 50 kpc. Bo should lie in the range [0.5 — 5] x 10~* 
gauss in order to produce UHECRs with the observed dis- 
tribution shown in this figure. 

It is worth emphasizing that this calculation was car- 
ried out without the use of several unknown factors often 
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Figure 7. Energy loss rate of a single electron in a turbulent time- 
varying magnetic field with Bq = 10~* gauss. In this diagram 
both the energy losses due to synchrotron and inverse Compton 
on the Radio and CMB photons in the acceleration region are 
shown. Unlike the case of protons and heavy ions, the radiation 
rate for an electron exceeds the acceleration rate in such a way 
that, for the given Bq, in a time of order of 10 years the electron 
will have lost all of its energy. 



required in approaches solving the hybrid Boltzmann equa- 
tion to obtain the phase-space distribution function for the 
particles. In addition, we remark that the acceleration mech- 
anism we have invoked here is sustained over 10 orders of 
magnitude in particle energy, beginning at 7 ~ 1; the UHE- 
CRs therefore emerge naturally from the physical condi- 
tions thought to be prevalent within the giant radio lobes 
of AGNs, without the introduction of any additional exotic 
physics (for a complete review of the bottom-up models, 
see Bhattacharjee & Sigl 2000 and other references cited 
therein) . 

To provide the possibility of observationally testing the 
model we have presented here, we show in Figure [7] the 
temporal evolution of the energy loss rate to first order in 
jhu/imeC^) (Blumenthal & Gould 1970) for a single elec- 
tron propagating through the same magnetic turbulence we 
have used to accelerate the protons and heavier ions. By 
estimating the flux of UHE protons Np escaping from one 
giant radio lobe, under the assumption of neutrality in the 
source, we can estimate the flux of accelerated electrons A^e- 
In principle, it is therefore possible to estimate the expected 
radio luminosity from these regions due to this particular 
acceleration process. 

We should point, however, that though a comparison 
of our results with the observations supports the viability 
of this model, our calculations are subject to several fac- 
tors we have not fully explored here. For example, the ob- 
served spectrum may b e affected by the cosmol ogical evolu- 
tion in source density (|Berezinskv et ahl (|2006l )). However, 
this omission will not be ove r ly con straining since a likeli- 
hood analysis l|Gelmini et al.l l(2007h ) of the dependence of 
the observed distribution on input parameters has already 
shown that, in the case of protons, for m ~ 0, where m is 
the evolution index in the source density, the HiRes obser- 
vations are compatible with a power-law injection spectrum 
with index —2.6. The anal ysis of the Auger data seems to 
confirm this (| Auger! l|2008d ) '). Thus, although source-density 



evolution may alter our results somewhat, our conclusions 
will probably not be greatly affected. In a more conserva- 
tive interpretation, the result presented here provides the 
injection spectrum from a single source. 

Second, we have not included the GZK effect for parti- 
cle energies above 50 EeV. This omission becomes progres- 
sively more important as the energy approaches lO'^" eV. 
These refinements, in addition to a more detailed analysis 
of the composition of primary UHECRs, will be reported in a 
forthcoming paper. Any discussion concerning the evolution 
of the instantaneous spatial diffusion coefficients parallel to 
and perpendicular to the background magnetic field, and on 
the transition to a diffusive regime, will also be reserved to 
a future publication. 

In our approach, we have also neglected the backreac- 
tion of the accelerated particles on the turbulent field which 
might increase the ratio |5f2|/|f2| , and bring about a possible 
local failure of the assumption of isotropy of the turbulence. 

We remark that the mechanism of stochastic accelera- 
tion presented here may be functioning even for a population 
of particles, protons or heavy nuclei, pre-accelerated to an 
initial energy E ~ 10^^ — 10^^ eV, e.g., by multi-shock fronts 
propagating at super- Alfvenic velocity. The corresponding 
gyroresonant wavenumber range in this case will decrease 

down to kmax/kmin ~ lO" — 10^. 



5 CONCLUSION 

We have shown that a region containing a Kolmogorov (tur- 
bulent) distribution of non-relativistic Alfven waves can ac- 
celerate particles to ultra-high energies. The physical pa- 
rameters in these regions are compatible with those believed 
to be operating in the radio lobes of AGNs. We have dis- 
cussed the predicted differential spectrum within the pa- 
rameter space of the model, characterized by the size TZ of 
the acceleration region and the turbulent magnetic energy. 
Possible tests of this model involve the synchrotron or IC 
emission by a population of similarly accelerated electrons. 

As the Auger observatory continues to gather more 
data, improving on the statistics, our UHECR source iden- 
tification will continue to get better. Eventually, we should 
be able to tell how significant the GZK effect really is, and 
whether the cutoff in the CR distribution is indeed due to 
propagation effects, or whether it is primarily the result of 
limitations in the acceleration itself. Given the fact that en- 
ergies as high as ~ eV may be reached within typical 
radio lobes, it is possible that both of these factors must be 
considered in future refinements of this work. 
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